function [CE,type_star,D_level_star,Pu_star,Pd_star] = ...
    Get_CE_RILA_optim_fun(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0)

%% Set model specifications & parameters
if Pu_Pd_same
    Pd_choices = 1;     % we will separately set Pd = Pu during the optimization.
end
F_num = length(F_choices);  B_num = length(B_choices);  Pu_num = length(Pu_choices);    Pd_num = length(Pd_choices);


%% State space
A_num = length(A_grid); A_grid_max = max(A_grid);   % grid of potential QRP account values A_t.


%% Matrizes to store results from dynamic optimization
EU_at_F = u(0.001) * ones(A_num,T,Pu_num,Pd_num);           % to store expected utility of being in current state (and time)
EU_at_B = u(0.001) * ones(A_num,T,Pu_num,Pd_num);           % to store expected utility of being in current state (and time)
EU_0_F = u(0.001) * ones(F_num,Pu_num,Pd_num);
EU_0_B = u(0.001) * ones(B_num,Pu_num,Pd_num);


%% Terminal utility
EU_at_F(:,T,:,:) = u(max(A_grid,0.01));
EU_at_B(:,T,:,:) = u(max(A_grid,0.01));


%% Recursion
if Floor_control
    for f = 1:F_num
        F = F_choices(f);
        for pu = 1:Pu_num
            Pu = Pu_choices(pu);
            for pd = 1:Pd_num
                if Pu_Pd_same
                    Pd = Pu;
                else
                    Pd = Pd_choices(pd);
                end

                Cap = FindCapRate_GHQ("Floor",F,Pu,Pd,r,tau,RILA_fee_rate,c_bar,R_S_GHQ_RN,weight);
                R_A = max(Pd*R_S_GHQ , -F) .* (R_S_GHQ<0) + min(Pu*R_S_GHQ,Cap) .* (R_S_GHQ>0);     % GHQ nodes for return credited to QRP account.

                for t = (T-1):-1:1
                    EU_tplus1 = EU_at_F(:,t+1,pu,pd);     % value function at end of period.
                    I_t = I_0 * (1+g_I)^t;  % PH's time-t contribution to QRP (no contribution at retirement time T).
                    for a = 1:A_num
                        A_t = A_grid(a);

                        A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max ); % GHQ nodes for QRP account value at year-end.
                        U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                          % GHQ nodes for utility at year-end.
                        EU_at_F(a,t,pu,pd) = weight' * U_tplus1;                                % expected utility of account values at end of period.
                    end
                end

                % Time 0
                EU_tplus1 = EU_at_F(:,1,pu,pd);     % value function at end of period.
                I_t = I_0;                          % PH's time-t contribution to QRP (no contribution at retirement time T).
                A_t = 0;
                A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max ); % GHQ nodes for QRP account value at year-end.
                U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                          % GHQ nodes for utility at year-end.
                EU_0_F(f,pu,pd) = weight' * U_tplus1;
            end
        end
    end
end
EU_star_0_F = max(max(max(EU_0_F)));

if Buffer_control
    for b = 1:B_num
        B = B_choices(b);
        for pu = 1:Pu_num
            Pu = Pu_choices(pu);
            for pd = 1:Pd_num
                if Pu_Pd_same
                    Pd = Pu;
                else
                    Pd = Pd_choices(pd);
                end

                Cap = FindCapRate_GHQ("Buffer",B,Pu,Pd,r,tau,RILA_fee_rate,c_bar,R_S_GHQ_RN,weight);
                R_A = min(Pd*R_S_GHQ+B , 0) .* (R_S_GHQ<0) + min(Pu*R_S_GHQ,Cap) .* (R_S_GHQ>0);     % GHQ nodes for return credited to QRP account.

                for t = (T-1):-1:1
                    EU_tplus1 = EU_at_B(:,t+1,pu,pd);     % value function at end of period.
                    I_t = I_0 * (1+g_I)^t;  % PH's time-t contribution to QRP (no contribution at retirement time T).
                    for a = 1:A_num
                        A_t = A_grid(a);

                        A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max );                % GHQ nodes for QRP account value at year-end.
                        U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                      % GHQ nodes for utility at year-end.
                        EU_at_B(a,t,pu,pd) = weight' * U_tplus1;                                 % expected utility of account values at end of period.
                    end
                end

                % Time 0
                EU_tplus1 = EU_at_B(:,1,pu,pd);     % value function at end of period.
                I_t = I_0;                          % PH's time-t contribution to QRP (no contribution at retirement time T).
                A_t = 0;
                A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max ); % GHQ nodes for QRP account value at year-end.
                U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                          % GHQ nodes for utility at year-end.
                EU_0_B(b,pu,pd) = weight' * U_tplus1;
            end
        end
    end
end
EU_star_0_B = max(max(max(EU_0_B)));


%% find combination of (F,B,Pu,Pd) that maximizes expected utility at time 0:
options = optimset('Display','none');
if EU_star_0_B > EU_star_0_F     % best Buffer is preferable to best Floor contract
    type_star = "Buffer";
    [b,pu,pd] = find(EU_0_B == EU_star_0_B);
    D_level_star = B_choices(b(1));
    EU_star_0 = EU_star_0_B;
else        % best Floor is preferable to best Buffer contract
    type_star = "Floor";
    [f,pu,pd] = find(EU_0_F == EU_star_0_F);
    D_level_star = F_choices(f(1));
    EU_star_0 = EU_star_0_F;
end
Pu_star = Pu_choices(pu(1));
Pd_star = Pd_choices(pd(1));

% Convert utility to Certainty Equivalent (CE):
CE = fsolve(@(ce)  1e+3 * (u(ce)/EU_star_0-1), 200,options);    % certainty equivalent to potential values of A_T.
